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Abstract 

The Allan variance was introduced fifty years ago for analyzing the stability of frequency standards. Beside its 
metrological interest, it may be also considered as an estimator of the large trends of the power spectral density (PSD) 
of frequency deviation. For instance, the Allan variance is able to discriminate different types of noise characterized 
by different power laws in the PSD. 

The Allan variance was also used in other fields than time and frequency metrology: for more than 20 years, it 
has been used in accelerometry, geophysics, geodesy, astrophysics and even...finances! However, it seems that up to 
now, it has been exclusively applied for time series analysis. We propose here to use the Allan variance onto spatial 
data. 

Interferometric synthetic aperture radar (InSAR) is used in geophysics to image ground displacements in space 
(over the SAR image spatial coverage) and in time thank to the regular SAR image acquisitions by dedicated satellites. 
The main limitation of the technique is the atmospheric disturbances that affect the radar signal while traveling from 
the sensor to the ground and back. In this paper, we propose to use the Allan variance for analyzing spatial data 
from InSAR measurements. The Allan variance was computed in XY mode as well as in radial mode for detecting 
different types of behavior for different space-scales, in the same way as the different types of noise versus the 
integration time in the classical time and frequency application. We found that radial AVAR is the more appropriate 
way to have an estimator insensitive to the spatial axis and we applied it on SAR data acquired over eastern Turkey 
for the period 2003-2011. Space AVAR allowed to well characterize noise features, classically found in InSAR such 
as phase decorrelation producing white noise or atmospheric delays, behaving like a random walk signal. We finally 
applied the space AVAR to an InSAR time series to detect when the geophysical signal, here the ground motion, 
emerges from the noise. 
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1. Introduction 

The Allan variance was designed by David Allan 50 years ago for characterizing the time stability of atomic 
clocks Q. This statistical tool quickly became a standard for researchers, engineers and manufacturers in the time 
and frequency field, to the extent that the Allan variance is not only considered as a stability estimator but as the 
stability itself! The use of this tool then spread to related fields such as inertial sensors Q, temperature 0 and 
voltage 0. Since 2000’s, the Allan variance has been used in many other domain including geodesy, to measure 
absolute displacement of geodetic stations 0, astrophysics, for selecting stable extragalactic compact radio sources 
from the VLBI program 0, geophysics, by monitoring the elevation of a volcano station to detect a change in 
the Allan deviation signature as a symptom of an eruption threat |[7|. The Allan variance is even used in finance 
for optimizing the financial returns 0 ! However, all these applications of the Allan variance deal with time series 
analysis. To our knowledge, this statistical estimator was never used over two-dimensional spatial data. This is 
the aim of this paper: modifying the 1-D Allan variance to a 2-D space Allan variance for fulfilling the needs of 
geophysics in terms of Earth surface deformation monitoring. 

Earth surface deformation is one of the most important observation to infer earth interior properties 
subduction zone activity 113 or natural hazard, including seismic hazard, landslide (B) or land subsidence 0 - 
While looking at large wavelength phenomena, measuring accurately how Earth deforms is a challenging task. A 
very effective way has been developed with geodetic techniques such as Global Navigation Satellite System (GNSS) 
or radar interferometry (InSAR). In both cases, measurements use the travel time of electromagnetic waves between 
satellites and Earth’s surface to derive the ground displacements. As Earth’s atmosphere is not neutral, taking into 
account the wave travel in this media is a key parameter to get reliable measurement of surface deformation. 

Here, we investigate the noise signature in InSAR and in particular the noise induced by the wave propagation into 
the atmosphere. InSAR has proven to be a very powerful tool to measure Earth surface displacements 0 Despite 
it measures the deformation in 1 dimension, this is the only technique offering large spatial coverage (hundreds 
of kilometers) and excellent spatial resolution (tens of meters). The main limitation is the atmospheric delays 
during the wave propagation in the atmosphere. Indeed, spatiotemporal variations of the air refractivity induce 
additional signal on the interferograms, that we called atmospheric phase screen (APS). Eollowing atmospheric 
conditions during the SAR images acquisitions, APS on interferograms can reach up to 10 cm, or more in case of 
extreme meteorological events. It is thus potentially much larger than the geophysical signal. A method to model 
the structure and, eventually, correct the APS is then an important goal to fully exploit InSAR as an accurate 
method for displacement measurement. The extension of AVAR to 2-dimensional images could be a powerful tool 
for achieving this purpose. 


H. Choice of the variance 

A. Expected properties 

The Allan variance is designed for assessing the large trends of the spectrum, not for detecting fine spectral lines. 
On the other hand, it can easily distinguish different power-laws in the spectrum of a time signal. This is exactly 
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Fig. 1. Example of white noise (left) and random walk (right). 


what we need in the analysis of InSAR images. 

We expect that the noise may be modeled as a power law of the wave number k, where k is the dual variable of 
the radius in polar coordinate, exactly as the frequency / is the dual variable of the time t in time series. Like in 
time and frequency, we could encounter a Gaussian white noise, of electronic origin or due to phase decorrelation, 
and an atmospheric noise exhibiting large scale correlation, likely with a power spectrum decreasing sls k~‘^ (random 
walk). Figure displays these two types of noise. Ideally, the response of Allan variance should be a plane with 
different slopes, depending on the type of noise: horizontal for the atmospheric noise and decreasing with the scale 
factor for white noise. 

This behavior should permit to separate the noise signature from the signal signature which depends on the 
topography of the images. The Allan variance must detect the large structure of the image signal. We ought then to 
identify the response of space AVAF[^for various features such as linear slope, a 1-dimensional sine, a 2-dimensional 
sine, with different orientations. Such patterns may be seen on Figure]^ 

Finally, on a time series of images covering the same area site over several years, we should be able to see the 
signal emerging from the noise as the date of the image increases. 

B. XY-AVAR 

The first attempt for extending the classical Allan variance from ID to 2D was the XF-AVAR. This variance 
relies upon a duplication of the double-stepped shape of AVAR along the first dimension, i.e. X-axis, in the second 
dimension, i.e. V-axis (see Figure |^. 

1) Definition: 

^ Since the classical Allan variance is often denoted AVAR, we will refer to the 2-diniensionaI space Allan variance as space AVAR. 
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Fig. 2. Example of an 1-dimensional sine image (left) and of 2-dimensional sine image (right). 




Fig. 3. From classical ID Allan variance to Xy-AVAR. 


a) Kernel: If XF-AVAR is considered as a 2D-filter, the kernel is its impulse response hxvix^y) which is 
the product of the kernel of classical AVAR along X by the kernel of classical AVAR along Y. Let us denote 
hx{x) the kernel of AVAR along X\ 

( 1 

hx{x) = —— if —Ax < X < 0 
v2Ax 

^ — if 0 < X <+Ax (i) 

v2Ax 

^ hx{x) = 0 if X < —Ax or x > +Ax 

where Ax is the scale factor along X. 

Similarly, let us denote hy (y) the kernel of AVAR along Y : 

/iy(2/) = -yX if -Ay<2/<0 

^ hviy) = if 0<y < +Ay (2) 

= 0 if X < —Xy or y > +Ay 


where Ay is the scale factor along Y. 
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Fig. 4. Squared modulus of the transfer function of XY-AVAR for Ax = 4 pixels and Xy = 3 pixels. 


The kernel of XY -AVAR is then: 


hxY{x,y) = hx{x) ■ hyiv). 


(3) 


b) Transfer function: The transfer function is the Fourier transform of the kernel (i.e. the impulse response). 
Moreover, since the kernel is the product of one function depending only on x by a function depending only on y, 
its Fourier transform is the product of the Fourier transform of these functions (see p. 139): 

Hxy{u,v) = Hx{u) • Hy{v) (4) 

where u and v are spatial frequencies in X and Y directions. 

Since the squared modulus of the transfer function of the classical Allan variance is given by (Tt): 


\HAVAR{f)\^ = 2 - 




(5) 


( 6 ) 


{xt/Y 

the squared modulus of the transfer function of XY -AVAR is then: 

\u ! m 2 .sin'^(7rAxw) sin^(7rAFf) 

c) First feeling: Figure shows the shape of this transfer function. This graph exhibits a grid of zeros since 
each time one of the functions Hx{u) or Hy{v) is equal to 0, it yields a line of 0 along the perpendicular axis. 
Therefore, XV-AVAR is strongly dependent on the axis orientation. This is clearly a disadvantage because the 
response of XV-AVAR may vary significantly versus the rotation of the image, although the image content is 
obviously the same. 

2 ) Responses of XY-AVAR: In order to test XY -AVAR, we applied it to the various noises and geometrical 


patterns defined in [ll-A 
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Fig. 5. Xy-AVAR of a 1-D sine (left) and of a 2-D sine (right). 


a) Response to a 1-D sine: The response of XF-AVAR to 1-D sine depends drastically on its orientation: if 
the sine is exactly in the direction of one axis, the XF-AVAR response is null; but this is no longer the case when 
the 1-D sine orientation is arbitrary as shown on Figure (left). 

b) Response to a 2-D sine: Similarly, if the 2-D sine are in the direction of the axis, XF-AVAR yields a null 
response. Otherwise, we can see an example of XF-AVAR response on Figure (right). 

c) Response to a drift: In this case also, the response may be null if the drift is oriented in the direction of 
one of the axis. 

d) Response to a Gaussian white noise: As expected, the response of XF-AVAR for a white noise is a plan 
decreasing proportionally to the scale factor. Figure (left) shows this plane seen almost edge-on. 

e) Response to a random walk: On the contrary, the response of XF-AVAR for a random walk is not at all 
flat (see Figure right): its curvature spans over 2 decades. Moreover, the surface exhibits a protruding central 
diagonal with flanged edges. This feature is not suitable because it could interfere with the signal signature. 

f) Searching for a noise giving a fiat horizontal response: The only noise which could give a flat horizontal 

response by XF-AVAR must have a power spectrum which decreases as rather than k~‘^ as a classical 

random walk. Figure [7] shows such a noise and the response of XF-AVAR for this noise. Firstly, we can see that 
this response is flat but not horizontal: it increases as the scale factor. Secondly, the response is only flat if the axis 
of the image are aligned to the axis of XF -AVAR. Moreover, the image of this noise is not at all natural: the X 
and F axis are very marked as a raster image. 


3) Conclusion about XY-AVAR: As expected in section ^II-Blc XF-AVAR is not suitable because it is too 
much dependent on its axis orientation. This has two consequences: 


some responses (for a drift, for a 1-D or 2-D sine) may completely vanish depending on the axis orientation 
a random walk does not yield a flat and horizontal response but a complex feature which could be confused 
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Fig. 6. Xy-AVAR of a white noise (left) and of a random walk (right). 




Fig. 7. Example of n ‘^v ^ noise (above) and the response of XY-AVAR for this noise (below). 


with the response for the signal. 

For these reasons, another space Allan variance algorithm must be sought. 

C. Radial AVAR: Allan’s hat 

We ought then to find a variance which does not depend too much on the orientation of the axis. Thus, we chose 
to define a kernel with a radial symmetry. The principle is to rotate the double-stepped shape of AVAR around 
the Z-axis (see Figure [^. In order to preserve the same surface for the central positive plateau and the peripheral 
negative one, it is necessary to respect a condition on the radii. Let us denote ri, the radius of the central positive 
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Fig. 8. From classical ID Allan variance to Allan’s hat. 


Step, r 2 , the radius of the peripheral negative step, Si the surface of the central positive plateau, and S 2 , the surface 
of the peripheral negative plateau: 


Si = Trrl = S2 = TT (rl - rl) ^ r2 = V2ri. 


(7) 


This leads to a 2-D shape that resembles a hat which we named the '‘Allan’s haf\ 

1) Definition: 

a) Kernel: Considering the full radial symmetry of Allan’s hat, the kernel may be written in polar coordinates 
as: 


^rad 


if 0 < r < ri 


\/27rr^ 

^rad(A —T if ri <r < V2ri 


( 8 ) 


V^Trr^ 

^rad(A^)=0 if r>V2ri 

where r = x‘^ -\-y‘^ and 0 = arctan^/x. Note that the factor l/\/2 has been kept for the sake of homogeneity 

with the definition of classical AVAR. As defined by ([^, radial AVAR is fully insensitive to 0. However, we must 
be able to detect the variations versus the scale factor Ax along X-axis as well as the variations versus the scale 
factor Ay along V-axis. For this purpose, instead of r, we define p as: 


P = 




\^X 


\\y 


(9) 


The Allan’hat is now elliptical: the semi-major axis of the central plateau is max {Ax, Ay}, its semi-minor axis is 
min {Ax, Ay}, its surface is Si = TrAxAy, and its bounds are given by p < 1. We can now rewrite the kernel in 
Cartesian coordinates: 

/ 

1 


V^TtAx Ay 
1 


^rad(^,^) = + 

^’■"d(^’^) = “V27rAxAK 

^rad(^,^)=0 if 


if 0 < 
if 1 < 






< 2 


( 10 ) 


> 2 . 
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Fig. 9. Squared modulus of the transfer function of rad-AVAR for Ax = 8 pixels and Xy = 6 pixels. 


b) Transfer function: In optics, it is well known that the Fourier transform of a circular aperture of transparency 
1 and of radius 1 is 2Ji(k)/k (the Airy pattern) where Ji{k) is the Bessel function of first kind and order 1 and 
k the wave number (see for example pp. 395-396). Therefore, the squared modulus of the transfer function 
of radial AVAR is: 




{'kXxXyY 


Ji(fc) _ Ji(V2fc) 

k 


( 11 ) 


Let us define the frequency variables u and v in such a way that k = We can rewrite |ll| as: 

ji Ji (^f{x\v? + xiv‘^)^ 


|-l^rad(^? I — 


(ttXxXy)^ 


( 12 ) 


^X\v? + Xlv‘^ f{X\v? + Xlv‘^) 

The shape of this transfer function may be seen on Figure 

c) First feeling: Rad-AVAR seems to be a promising algorithm. Let us study its responses for the various 


noises and geometrical patterns defined in ^ II-A 
2) Responses: 

a) Response to a ID sine: The response to a ID-sine changes slightly with the axis orientation (basically, for 
a 90° rotation angle, the behavior along the A-scale factor is exchanged with the one along the V-scale factor and 
conversely) but the response is never identically null (see Figure p0| ). 

b) Response to a 2D sine: Similarly, the response to a 2D-sine changes slightly with the axis orientation but 
the response is never identically null (see Figure pT] ). 

c) Response to a linear drift: The response passes from a plane increasing proportionally to the A-scale factor 
when the slope of the drift is oriented along the A-axis to a plane increasing proportionally to the V-scale factor 
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Angle: 0° Angle: 45° 



Fig. 10. Response of Rad-AVAR to a ID sine with different orientations. 




Angle: AS'’ 




Fig. 11. Response of Rad-AVAR to a 2D sine with different orientations. 


when the slope is oriented along the F-axis (see Figure This is then an important property of Rad-AVAR. 

d) Response to a Gaussian white noise: The response of Rad-AVAR is a plane decreasing with the scale factor 
as for XY -AVAR. Here also, it is almost seen edge-on. 

e) Response to a Gaussian random walk: For a random walk, Rad-AVAR gives a quasi flat and horizontal 
plane: the Z-amplitude ranges from -0.9 to -0.7 (from 0.13 to 0.20 on a linear scale) whereas the X and Y- 
scale factors range over 2 decades. The residual fluctuations of the response are due to the difference between the 
mathematical expectation of Rad-AVAR, which is exactly flat and horizontal, and the estimator of this variance 
which remains a random variable (as for the classical Allan variance). Nevertheless, its fluctuations are very 
limited and this tool is then entirely able to distinguish between a white noise, a random walk and a geophysical 
signal. 




















































































Fig. 12. Response of Rad-AVAR to a drift with different orientations. 



Fig. 13. Response of Rad-AVAR to white noise (left) and to a random walk (right). 


3) Conclusion about radial AVAR: Radial AVAR has exactly the expected properties defined in ^II-A It is then 
perfectly suitable for our purpose. 


III. Application to InSAR data 

A. Data and method 

We worked with synthetic aperture radar (SAR) images acquired by Envisat, satellite from the European Space 
Agency (ESA). This dataset has been primarily used in 03 to map the surface deformation in eastern Anatolia. SAR 
images are acquired in illuminating the ground with electromagnetic waves (microwave frequency) and collecting 
the reflected signal. This remote sensing technique uses the round trip travel time and amplitude of the waves to 
determine physical properties of the targets on the ground that refiected the signal. Two quantities are measured : 
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(i) the amplitude of the reflected signal that corresponds to the energy sent back in the satellite direction and (ii) the 
phase of the reflected wave. Interferograms are generated in differentiating the phase of two SAR images acquired 
at different time pQ| . Considering that the ground response are identical for both SAR images (i.e, these terms 
cancel each other out in the phase substraction operation), the observed interferometric phase (p can be written as 
the sum of four terms (ID We add, here, one more term to distinguish explicitly between two kinds of atmospheric 
delays and the other sources of noise : 


^ — ^def ^orb ^strat ^turh ^noise (13) 

where cpdef is the phase component related to ground movements, (forh is the residual orbital phase, (ptropo 
is the stratifled atmospheric phase delay correlated with elevation, (fturh is the turbulent delays due to the 3- 
dimensional refractivity distribution of the atmosphere, and noise is the remaining term that includes instrument 
noise, DEM errors, or phase decorrelation. Phase decorrelation occurs when ground surface changes between two 
SAR acquisitions (vegetation growth, ground erosion) or by geometrical decorrelation when the perpendicular 
baseline between the two images is too large. When decorrelation happens, wrapped phase takes random values 
between —tt and tt, that can be assimilated as white noise. Figure shows phase decorrelation. We obtain a white 
noise space AVAR signature as obtained in our previous test (Figure [^. Such interferograms are then discarded 
for ground displacement measurements. 

In our data processing we corrected interferograms for stratifled tropospheric delay, strata using ERAI, a global 
atmospheric model p^ . This correction is a key step to obtaining a successful measurement of tectonic signals 
by InSAR, although this correction is not always successful and relies on the accuracy of the global atmospheric 
model. shows that ERAI does not predict properly the tropospheric delay for 10% of the case, ip orb is mostly 
caused by a drift in the local oscillator onboard the Envisat satellite p4| . This drift only causes a ramp in the 
across-track direction (range) of the interferograms and we therefore made no ramp corrections in the along-track 
direction (azimuth), to avoid removing tectonic signals related to the motion of the Anatolian and Arabian plates, we 
estimate the residual orbital signal by an east-west ramp across each interferogram only on their Eurasian part that 
is considered as stable. Others sources of noise can be turbulent atmospheric delays (fturb^ DEM errors or thermal 
noise, ipnoise- Note that, for the selected interferograms, (fturb is much larger than ipnoise- Figure p~5] shows single 
interferograms containing both geophysical signal and noise. As the ground movement is small, noise dominates 
the scene. Turbulent atmosphere patterns decorrelate within 24 hours | [25| . So, as SAR images are taken at different 
dates (shortest repetitive acquisition is 35 days), interferograms that do not share common images, will exhibit 
different noise signal patterns (Figure p3]). 

To derive ground displacement rates from noisy interferograms, used small baseline time series analysis 
This method uses interferograms, computed with the same image dataset, to retrieve phase increment between 
two images acquired successively in time. As said previously, most of the remaining interferogram noise is high 
frequency in time. A smoothening operator is then applied to limit phase variations due to turbulent atmospheric 
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Fig. 14. Right: Interferogram with phase decorrelation (red square). This occurs when spatial baseline is too large between both SAR images 
or when surface properties (mostly vegetation and erosion) changed between both SAR images acquisition. Left : Space AVAR of the phase 
values located in the red square. This response corresponds typically to a white noise. 


delays. This step allows to estimate the surface displacement history of each coherent pixel from the set of small- 
baseline interferograms. For tectonic application, p9| computed the linear component of the time series result for 
each pixel, to derive the steady surface displacement velocities. This result can be explained by simple elastic model 
where plate borders are represented by buried dislocations, creeping at a rate allowing to accommodate relative plate 
motions (i.e, here Eurasia, Anatolia and Arabia). Figure shows observed (by InSAR) and modeled velocities 
maps. Note that maps are shown here in the geographic referential and the phase has been converted in centimeters 
(for Envisat satellite, In = 2.8 cm). For simplicity reasons, we will work with radar coordinates in the following 
and phase amplitudes will be in radian. 

We tested the Allan Variance on various InSAR items. First on synthetic tests, to see how the 2D-AVAR reacts, 
and then on the real InSAR data presented above. 

B. InSAR signal modeling 

1) Tectonic signal modeling: To have an idea of how the 2D-AVAR of the geophysical signal evolves with 
time, we modeled the tectonic signal in radar coordinates for each SAR image date. The spatial AVAR is then 
computed at each time step. The first SAR acquisition of the time series is our reference image. There is thus no 
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Fig. 15. Examples of interferograms obtained along track 493. Colorbars indicate phase changes in radians, a) Typical atmospheric patterns due 
to the turbulent and spatial heterogeneities of the air refractivity. b) On top of the atmospheric delays, higher frequency signal is superimposed 
due to phase decorrelation, c) Example of miscorrected interferogram where the effect of the stratified atmosphere is still visible (black circle), 
d) Interferogram with residual “orbital ramp” along the range axis. 


deformation. We will then consider the deformation between the date of each other images and this reference. We 
assume the ground displacement to be linear in time. So, knowing the plate velocities, it is then possible to compute 
the expected geophysical signal between each SAR image and the reference image (Figure [17^). As the tectonic 
signal has a large spatial wavelength content and increases linearly with time, we expect the AVAR to increase with 
time and as the scale factor increases. 

Tectonic signal is mostly a large wavelength signal. Applying the spatial AVAR on the modeled signal, we see 
that the variance increases both in azimuth and range axis. However, du to the geometry and movement of the 
plates, AVAR reaches larger values along the azimuth axis (Figure [TT])). As the signal grows in time, the shape and 
AVAR range do not change, but the maximum AVAR value, found for large scale factor, increases (Figure [T^). 

2) Atmospheric noise modeling: Based on Tatarski’s work | [27| , notably, Hanssen proposes a new atmospheric 
model for radar interferometry considering three regimes of atmospheric delays as observed by InSAR 1^ . Each 
of these regimes corresponds to a particular scale factor and is characterized by an exponent of the power spectrum 

I scale-factor > 2 km: P(f){k) oc 

II 0.5 < scale-factor < 2 km: P(f){k) oc k~^/^ 

























IEEE TRANSACTIONS ON ULTRASONICS, EERROELECTRICS, AND EREQUENCY CONTROL, 201 


16 



Fig. 16. Observed and modeled ground velocities, (a) Observed horizontal velocities from InSAR time series analysis in mm/yr parallel to 
the satellite line of sight (LOS) direction (black and white arrows show the LOS and flight directions, respectively) and in an Eurasia-flxed 
reference frame. Velocities are shown positive in the east-southeast direction (103 N). (b) Elastic back slip model modeling. The black cross in 
(a) points the location of the time series of Figure 


III 0.01 < scale-factor < 0.5 km: P(f){k) oc k 

where k is the wave number and P(f){k) the power spectrum. This model assumes isotropy and no ionospheric 
influences, that is a verifled for C-band SAR data as used in this study. However, the model does not take into 
account topography related delays, (ftropo^ what might be a problem if this term is not properly removed. 

We can, thus, model the atmospheric noise by the following dependence on the wave number: 

where 

• ki 2 is the cut-off wave number between regimes I and II: ki 2 = 0.5 km“^ 

• is a free parameter denoting the general level of atmospheric noise 

• /13 is the level of the regime III which value is flxed by the previous parameter: since the limit between regimes 
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Fig. 17. (a) Modelled tectonic signal in radar coordinates corresponding to the line of sight (LOS) ground displacement, (b) Space Allan 

variance computed from (a), (c) Evolution of logiO(AVAR) maximum value for each SAR acquisition date. 


II and III corresponds to scale-factor of 0.5 km, i.e. /C 23 = 2 km“^, it yields 

— ^3^23^^ ^ hs = hk2^ = hj 4. 

The atmospheric noise model is then: 




hk + A:-2/3/4) 


(14) 


A:+ 0.5 

where k is expressed in km“^. 

In our case, a pixel corresponds to 640 m and an image contains typically 600 pixels. Thus, k ranges between 
0.0025 and 0.8 km“^ (Nyquist frequency). 

shows the Hanssen’s model. It is very interesting to notice that this model remains close to a k~‘^ 


Figure 


18 


asymptote within the A:-range which corresponds to our images. It means that the response of the atmospheric noise 
to space AVAR should be very similar to a plan almost flat and horizontal which confirms our intuition of §II-A 


This model may be used for simulating images of atmospheric noise. Such an image is displayed in Figure 


19 The space Allan variance of this sort of images may thus be computed. Figure ^ exhibits the average of 


the response of AVAR for 25 uncorrelated simulated atmospheric noises like Figure Averaging 25 responses 
cancels out the deviation of each individual response from the theoretical variance of space AVAR. Thus, Figure 


20 is representative of the response of space AVAR for an atmospheric noise following the Hanssen’s model. 


As expected, taking into account the very thin range of Figure 2 ;-scale, we can see that the response of space 
AVAR is almost horizontal and fiat. More precisely, it looks like a cupola: horizontal near the center and decreasing 


on the edges. For going further, we must analyze Figure 21 which shows the power spectrum of the atmospheric 
noise multiplied by k‘^, in order to identify the areas where the slope follows a A:“^ asymptote rendered as an 
horizontal curve. Indeed, this figure shows that the atmospheric power spectrum follows exactly 3 . k~^ asymptote 
around A: = 0.3 km“^, yielding a fiat response of space AVAR for scale factors around 5 (0.7 on log 10-scale). For 















P^k) / h 
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Wave number k / km'^ 


Fig. 18. Model of atmospheric noise. The three regimes are visible: regime I at left (for k < 0.05 km“^, dashed red asymptote), regime 
II near the center (k ^ 0.5 km“^, dashed green asymptote), regime III at right (for k > 5 km“^, dashed purple asymptote). However, the 
atmospheric noise model remains relatively close to the k~‘^ asymptote (dash-dotted cyan line). 
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Fig. 19. Simulation of an atmospheric noise according to the Hanssen’s model described by relationship (14}. 


larger k, the curve of Figure exhibits a decreasing slope: this means that the power-law of the power spectrum 
is lower than —2 yielding a positive slope for space AVAR at small scale-factors. On the other hand, for k smaller 
than 0.3, the curve of Figure [21] exhibits an increasing slope: the power-law of the power spectrum is then higher 
than —2 and space AVAR decreases as the scale factor increases. 

We can now compare this theoretical model with single interferograms acquired over the area presented on 


Figure For that purpose, we chose interferograms with small temporal baseline (i.e small timespan between 
the two SAR acquisitions) in order that the tectonic signal is negligible with respect to the atmospheric delays. 
Space AVAR would thus reflect the atmospheric structure. Figure shows that typical space AVAR pattern for 

a corrected interferogram (i.e (fstrat and ^orbit have been removed) differs from Hanssen’s model at large scale 
factors. Space AVAR response of the digital elevation model (DEM) also displays a quasi-plane increasing with 
scale factors. So even, if DEM related components have been removed, residual of this component might still be 
strong enough to contaminate the space AVAR of the APS response. Space AVAR shows thus that in areas where 
strong elevation variations take place, like in Turkey, the influence of the topography cannot be neglected to model 
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Fig. 20. Response of space AVAR for an atmospheric noise. This response is the average of the response of AVAR for 25 uncorrelated simulated 
atmospheric noises (Figure [T9| is on these). It must be noticed that z-scale ranges from to i.e. from 3 • 10“^ to 7 • 10“^, 

which proves that this response is nearly flat. 


the atmosphere for InSAR. However, this effect must be relativized since the magnitude of this increasing of space 
AVAR for high scale factors is limited to half a decade (from to 10“^, i.e from 3-10“^ to 10“^). 

An empirical way has also been proposed to model the atmospheric topology for InSAR p9| . This method 
estimates directly the variances and autocovariances of interferograms. For this purpose, semi-variograms j(h) is 
used to estimate the InSAR variances and covariograms C(h) is used to estimate the spatial correlation in the data 
| [30| . They are defined as followed : 


i,3 

d(i,j) = h 


(15) 


where ip{i) and ip{j) are the interferometric phases of pixels i and j, respectively, is the distance beween 

i and j, and N is the number of pixel pairs such that d{i,j) = r. Similarly, one can define the covariogram as: 





P^kj.k" / h 
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Fig. 21. Power spectrum of the atmospheric noise model multiplied by k‘^ within the observable range. The k‘^ factor allows us to clearly 
identify that the power spectrum follows 3l k~‘^ power-law (horizontal asymptote) around k = 0.3, a power-law with an exponent > —2 for 
lower k, and a power-law with an exponent < — 2 for higher k. 


■ ‘^ 0 ) ( 16 ) 

i,3 

d{i,j) = h 

From the semi-variograms, one can estimate the variance of the interferograms when the function flattens (Figure 
[22| ). Covariogram indicates the correlation length in the data and how fast it decorrelates. In case of white noise 
(i.e no correlation), covariogram function would be constant-zero function and semi-variogram would be a constant 
function at the variance value. Covariagram function can be estimated by a an exponential decay b • exp(—/i/a). 
The exponential decay function is then used to generate random, similar spectral, noise (Figure [^ . 

C. Real data 

We now applied the 2D-AVAR on data acquired from Envisat satellite over an area covering the eastern part of 
the Anatolian plate. We combined 18 SAR images into 32 interferograms. Figure shows image acquisitions date 
as a function of their relative position (perpendicular baseline). Note that to keep a good phase coherence over the 
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Fig. 22. (a) Interferogram computed with images acquired with a 35-day timespan. (b) Semi-variogram, ^{h) and (c) covariogram, C{h), 

computed for this interferogram. 


interferogram, perpendicular and temporal baseline between two images is limited. Based on this criteria, we find 
the best network that connect each image. Each interferogram is corrected for ip orb and Pstrat p9| - We computed 
for each interferogram the spatial AVAR before and after the correction. We observe that the correction step clearly 
decreases the mean AVAR, confirming that large wavelength have been reduced. The next step is to inverse the set 
of interferograms to establish a time series of the phase evolution. 

D. Data Inversion 

1) Inversion without smoothing: Corrected interferograms are then inverted to solve for the incremental de¬ 
formation between two successive images using a least squares inversion method (D For each pixel, treated 
independently from its neighbours, we solve: 

d = G.m (17) 

where d includes N interferograms observations, m corresponds to M-1 incremental displacements between the M 
time steps, and G is a by M — 1 matrix containing zeros and ones, based on the stating that the interferometric 
phase, pij, is the sum of successive phase increments between image i and image j : pij = Images are 

ordered by dates. 


The inversion uses the redundant information of the 32 interferograms to reduce some errors due to interferogram 
formation. 

2) Inversion with smoothing: To decrease atmospheric delay or APS, noted pnoise m equationwe introduce 


temporal smoothing in the inversion as another constraint, by minimizing the curvature of the inverted phase temporal 
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Fig. 23. (a) Same interferogram as in Figure [2^ (b) corresponding space AVAR diagram, (c) Modeled atmospheric noise based on the 

interferogram covariogram and (d) corresponding space AVAR diagram. 


evolution 0- The system becomes : 



d, G, and b are weighted by the matrix W (d' = Wd, G' = WG, and b' = Wb). 7 is the smoothing coefficient 
introduced to ponderate the minimum curvature constraints ^ (where t is time). ^ | • is evaluated with a 5 points 
finite difference scheme centered on each acquisition date i. The weight, W, applied on each interferogram is the 
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Fig. 24. Diagram of Envisat images for descending track 493. Relative perpendicular baselines are plotted as a function of acquisition dates. 
The 18 SAR images are combined into 32 interferograms (grey segments). 


product of two terms : (1) the first term equalizes the weight of all images in the inversion, (2) the second term 
characterizes the interferograms “quality”, qij, defined as the inverse sum of both images i and j APS amplitudes. 


3) Times series results and spatial AVAR: Figure [25] shows the time series of the ground displacement in the 
satellite line of sight (LOS) with or without smoothing in the inversion. We see that after smoothing, ground 
displacement rate becomes steady. However, it does not mean that all the high frequency APS have been removed. 
Actually, if it seems so in time for one pixel, the resulting smoothed displacement rate may change from one pixel 
to another due to different APS conditions. Spatially, APS are thus not fully removed, but its amplitude clearly 
decreased. 

As the tectonic signal contains long wavelength, to characterize it, we took the maximum value of the AVAR, 
maxav obtained at largest scale factor. We did that for the non-smoothed and smoothed times series. Figure ^6 


shows the evolution maXav with time (as a reminder, in times series, APS component amplitude is decorrelated 
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Fig. 25. Example of ground displacement times series for a pixel located on the Anatolian plate (see the cross in Figure for location). 
Black and red lines represent the ground displacement between 2003 and 2011 without or with smoothing, respectively. 


from one to another SAR image, and depends on specific atmospheric conditions during the image acquisition, 
while tectonic signal grows steadily with time). We see that when no smoothing is applied, no particular behavior 
of maxav shows up with time. It seems then that APS dominates the time series signal. On the other hand, for 
the smoothed time series, maXav increases linearly with time after 2 years. It seems then that the geological signal 
emerges from the noise after 2 years. Earlier, plate motion signal is in competition with the APS. Figure [27] shows 
the time series signal after 1 year, where the small geophysical signal is lost in the noise, and the final signal over 
the ^7-year period of observation. 
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Time from the reference image [yr] Time from the reference image [yr] 

Fig. 26. Maximum of logio(AVAR) value found at each incremental time step for ground displacement time series computed (a) without 
smoothing or (b) with smoothing. On (a), APS clearly dominates the logioAVAR maximum value. On the contrary, on (b), after 2 years, when 
the tectonic signal starts to be strong enough, logioAVAR maximum values increases proportionally with time, as for the model (Figure [TtJ:). 


E. Conclusion about InSAR application 

Space AVAR allows to distinguish between classical signals that one can deal with : white noise, random walk 
and others. In InSAR, different signals sum up and it can be difficult to observe the geophysical signal, typically 
here ground surface displacement due to plate tectonics, among noise sources such as atmospheric delays or phase 
decorrelation. We saw that the phase decorrelation is easily identifiable and responds as white noise. Based on 
theoretical, model atmospheric signature is also characteristic and differs significantly from other signals. However, 
in analyzing our data, we noted that atmospheric delays had additional large wavelength signal. This could be due 
to large topographic variations that would induce a correlated elevation component, not taking into account in the 
atmospheric model. Moreover, knowing the geophysical signal, we were able to compute its space AVAR signature 
which one feature is the larger AVAR values at large scale factors. As the tectonic signal is linear in time, it is 
easily verifiable if the time series has this characteristic. However, it is more difficult to assess the geophysical 
signal in space when it is mixed with different type of noise. In taking the space AVAR value for large scale factor, 
we were then able to identify when the geophysical signal emerged distinctly from the noise. 

IV. Conclusion 

This study shows an original way to compute and use the Allan variance in the space domain. It turns out 
that radial AVAR, called here Allan’s hat, is the best way to obtain a stable AVAR solution, insensitive to the axis 
orientation. A fast algorithm has been elaborated in computing the Fourier transform transfer function of the Allan’s 
hat, that corresponds to a sum of two Bessel functions divided by the radial wave number (a sum of two Airy 
patterns). 
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We then applied this formalism to compute the space AVAR for InSAR data. It allowed to retrieve the main 
characteristics of noises that contaminate the ground displacement measurements obtained by computing the phase 
difference of two SAR images. We could clearly observe white noise due to phase decorrelation or quasi-planar 
AVAR features due to atmospheric delays. To distinguish between noises and geophysical signal, we extracted the 
AVAR value for large scale factors, where the geophysical signal is the strongest. It showed great stability and 
revealed the ongoing geophysical signal that emerged from noise after a period where noise and tectonic signal 
were of comparable amplitude. 

Through InSAR data, we showed that space AVAR can be useful to detect noise and characterize 2-D signals. 
It can thus be used to a multiple of other domains that deal with spatial data. Compared to classical AVAR, space 
AVAR may give more complicated features to analyze, especially when occasionally the AVAR surface is bumpy 
and not planar, but seems very efficient for distinguishing noises from signal. Space AVAR could certainly be used 
more generally for image processing in other fields or even for other 2-D signals (surface analysis,...). Moreover, 
the technique used here to go from 1-D to 2-D can also be extended to 3-D. We can imagine applications like 
volume temperature monitoring in a room, 3-D atmospheric evolution, etc. 
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Fig. 27. (a) Signal obtained after 1.05 years of observation and (b) the corresponding space AVAR diagram, (c) Signal obtained over the 

whole time series and (d) the the corresponding space AVAR diagram. 























